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ABSTRACT 

We demonstrate that the sensitivity of high-precision pulsar timing experiments 
will be ultimately limited by the broadband intensity modulation that is intrinsic 
to the pulsar's stochastic radio signal. That is, as the peak flux of the pulsar ap- 
proaches that of the system equivalent flux density, neither greater antenna gain nor 
increased instrumental bandwidth will improve timing precision. These conclusions 
proceed from an analysis of the covariance matrix used to characterise residual pulse 
profile fluctuations following the template matching procedure for arrival time estima- 
tion. We perform such an analysis on 25 hours of high-precision timing observations 
of the closest and brightest millisecond pulsar, PSR J0437— 4715. In these data, the 
standard deviation of the post-fit arrival time residuals is approximately four times 
greater than that predicted by considering the system equivalent flux density, mean 
pulsar flux and the effective width of the pulsed emission. We develop a technique 
based on principal component analysis to mitigate the effects of shape variations on 
arrival time estimation and demonstrate its validity using a number of illustrative 
simulations. When applied to our observations, the method reduces arrival time resid- 
ual noise by approximately 20%. We conclude that, owing primarily to the intrinsic 
variability of the radio emission from PSR J0437— 4715 at 20 cm, timing precision in 
this observing band better than 30 - 40 ns in one hour is highly unlikely, regardless 
of future improvements in antenna gain or instrumental bandwidth. We describe the 
intrinsic variability of the pulsar signal as stochastic wideband impulse modulated 
self- noise (SWIMS) and argue that SWIMS will likely limit the timing precision of 
every millisecond pulsar currently observed by Pulsar Timing Array projects as larger 
and more sensitive antennas are built in the coming decades. 
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1 INTRODUCTION 

The most fundamental property of radio pulsars is their 
periodic series of radio pulses that enable their discovery 
and a myriad of timing applications. A sub-class of pul- 
sars, known as the millisecond and recycled pulsars, have 
spin periods between 1.4 and a few tens of ms and typical 
spin-down rates of P ~ 10~'^°. Their short periods and low 
braking torques make them especially good c locks, and these 
pulsa rs exhibit the highest timing precision jMatsakis et al.l 
1 19971 ). For most of these pulsars a simple model of the pulsar 
spin-down, astrometric and orbital parameters can be deter- 
mined, enabling the mean time-of-arrival (ToA) of pulses to 
be predicted accurately and precisely. These can be used 
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for several applicat ions, such as testing the general the- 
ory of relativity fe.g. lTavlor fc Weisbergiri982l : lKramer et al.l 
L2006b|) , detecting irregularities in terestrial time stan dards 
CPetit_ & TavcUa 1996; Rodin 2008; Hobbs et al.ll2010l ') and 
to attempt the first direct detection of a stochastic back- 



ground of gravitational waves (see, e.s 


j^. HellinES & Downs 


19831: iFoster & Backed 1990l: DemorestI 


boOTi: lYardlev et al. 


2010l:lvan Haeisteren et al. 201ll). The wealth of information 



already derived from the precision timing of millisecond ra- 
dio pulsars has led many authors to predict the kind of pul- 
sar timing science possible with the Square Kilometre Array 
(SKA) by linearly extrapolating current telescope sensitivi- 
ties to that of the SKA. 



The closest a nd brightest m i lliseco nd pulsar known, 
PSR J0437-4715 l|johnston et all 119931 ) has been studied 
by numerous authors with steadily improving instrumen- 
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tation. ISandhu et alj l| 19971 ) observed the pulsar using an 
autocorrelation spectrometer with 128 MHz of bandwidth, 
and could model pulse arrival times over two years with a 
post-fit residual (the difference between observed and pre- 
dicted arrival times after fitting for the pulsar s pin, astro- 
metric and binary parameters, etc.; lTavloilll99j ) standard 
deviation of 500 ns. Noting that the formal uncertainty of 
arrival time estimates was typically around 50 ns, the au- 
thors concluded that the ir result s were limited by polari- 
metric calibration errors. iBrittonI l|2000l ) first proposed the 
use of the Stokes invariant interval to mitigate the prob- 
lems caused b y polarisation calibration . This was later im- 
plemented bv lvan Straten et al. I (2001), who used a combi- 
nation of 16 MHz and 20 MHz baseband recording systems, 
typical integrations of 1 hour duration, and coherent dedis- 
persion to obtain a root-mean-squaxe (rms) timing resid- 
ual of 130 ns over 3.4 yr. Us ing new and i mpro v ed ca libra- 
tion methods developed by Ivan StratenI (|2004 120061 ) and 
a new baseband recording and processing system (CPSR2; 
BailesI l2003l: iHotanI l2006f ) with 128 MHz of bandwidth, 
Verbiest et al.l ~ 200^ ) achieved 199 ns over 10 years. 

None of the above studies achieve the timing precision 
predicted by the formal uncertainty in arrival time esti- 
mates. When observing PSR J0437-4715 in the 20 cm band 
at the Parkes 64m observatory, the expected rms timing 
residual from a 256 MHz band with 21 K system temper- 
ature is about 10 and 80 ns in one hour and one minute 
of integration, respectively. These uncertainties are derived 
from the template- matching method used for pulsar timing, 
in which each observation of the average pulse profilsQ 0{t) 
is modelled as a scaled (A) and offset (B) template S{t), 
rotate d by some phase shift with additional white noise 
N{t) (|Tavloij|l992l : lBailesll201ol ). 



0{t) = AS{t - (l>) + B + N{t). 



(1) 



It is generally assumed that the summation of many hun- 
dreds or thousands of pulses leads t o a stable pulse pro - 
file that is characteristic of the pulsar (jHelfand et al.|[l975h . 
Consideration of only additional white noise, N{t), in the 
above equation is equivalent to assuming that the sys- 
tem equivalent flux density is the only signiflcant source 
of noise. However for bright sources and/or high gain an- 
tenna, this assumption is no longer tenable in at least two 
circumstances. Firstly, when the flux density of the pulsar 
approaches the system equivalent ffux density (SEFD) of 
the receiver additional noise propo rtional to the pulsar's 
ffux density beco mes significant (e .g. [Kulkarni" 1989*; 'Gwinn! 
2001, 2004, 2006; van Stratenll2009l :lGwinn & Johnson 2011; 
Gwinn et al.ll201li r Secondly, it is known that each single 



pulse can have very di fferent morphology and can occur 
at di f ferent pulse phase jPrake fc Craft|[l968l : iHelfand et al.l 
1 19751 : iJenet et al.|[l998l ). Even after integrating over many 
pulse periods, this subpulse modulatioio can introduce de- 
tectable ffuctuations in the average profile shape and thereby 
contribute additional noise to timing data. We discuss the 



^ throughout the paper we refer to the observed averaged phase 
resolved light curve of the pulsar as the pulse profile or pulse 

shape 

^ The term "phase jitter" (e.g. lCordes fc Shannonll201Cll ) is some- 
times used to describe this phenomenon but we find it somewhat 
misleading as it is not only the pulse phase that varies. 



noise balance in more detail in ij2l and argue that these two 
contributions should be considered together as they are re- 
lated and are described by the same statistical model. 

We note that the presence and importance of pulse 
profile variability have been discussed in many different 
contexts. Some "classical" pulsars have been observed to 
change between two or more stable profiles - a phenomenon 
known as mode changing - on time-scal es of minutes to 
hours Ce.g. lBackeilll97d iBai-tel et al.lll982l ). On longer time 
scales, pulsars have been discovered whose emission com- 
plete ly switches off for many days, weeks or e ven months 
(e.g. iDurdin et all 1 19791 : iKramer et al] l2006al ). Recently, 
iLvne et al.l (|2010l) have shown that the pulse profiles for 
many pulsars switch be tween two unique states on time- 
scales of months to vears. lKarastergiou et all l|201ll ) recently 
detected a transient component in PSR J0738— 4042, vary- 
ing on time-scale of years or decades. 

The connection between the pulse shape changes and 
timing noise was made soon after t he discovery of timin g 
noise in the pulsar observations by iBovnton et al] (| 19721 ). 
They studied optical timing observations of the Crab pulsar 
and discovered a noise component in the timing residuals 
which was well modelled as a random walk in the pulsar 
spin frequency. The authors also considered a random walk 
in the pulse phase and spin frequency derivative, but found 
no evidence of such noise in th eir data. Th is analysis was 
extended in a series of papers (|Grothl ri975a .h:.c) . The au- 
thor presented an analysis method suitable for studies of 
data with inherent timing noise. This improved methodol- 
ogy led to the conclusion that noise in the Crab pulsar tim- 
ing is dominated by a random walk in the spin frequency 
but a random walk in pulsar phase m i ght al so be present. In 
the meantime, [Manchester fc Tavlorl (| 19741 ) described tim- 
ing noise for two slow pulsars from radio observations. A 
few years later, another series of pape rs (Hclfand ct al. 198^; 
ICordej|l98AlCordes fc Helfandlll980l ) presented statistics of 
timing noise for 37 bright pulsars and concluded that it is a 
ubiquitous phenomenon. These authors presented a careful 
analytical description of random walks in pulsar phase, spin 
frequency and its derivative and the uncertainties in the es- 
timation of their parameters. The last paper in the series 
pointed out that the random walk in the pulsar phase can 
be due to the random pulse shape changes but concluded 
that it was unlikely to be the dominant source of timing 
noise in the available d ata. A different da t aset w as analysed 
in a similar manner bv lCordes fc Downj (|l985l ) who stated 
that either excessive jitter or pulse shape changes are im- 
portant for a significant fraction of their sample. They also 
pointed out that the pulse shape changes are likely to be 
univer sal but t heir im portance varies from object to object. 
Later, ICorded l| 19931 ) detected pulse shape variability in 11 
out of 14 studied objects. These variations were consistent 
with being ca used by summatio n of a finite number of pulses. 
A year later [Kastai et al.l (|l994l ) studied two millisecond ra- 
dio pulsars and discovered timing noise in one of them. The 
general continuity of properties between classical and mil- 
lisecond pulsars suggests that pulse shape changes may be 
common in millisecond pulsars as well. 

The profile variability of millisecond pulsars (MSPs) 
has been studied in relatively few cases. The single pulses 
from PSR J1939-I-2134 show no s ubpulse structure over 
selected ranges of pulse longitude l|jenet et al.l 1200 ll ) but 
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emit giant pulses as much as 300 times brighter than the 
average pulse, that are narrower and systematically de- 
layed with respect to the main and interp ulsc components 
l|Cognard et al.lll996l : iKinkhabwala fc Thorset t 2000). Sev- 
eral other groups have argued that MSPs exhibit pro- 
file shape changes. Some are associated with different 
viewing geometries or with gravitational spin-o r bit cou- 
pling; e.g. PSR B1913-H6 (IWeisberg ct al.' 1989: iKramej 



1998 ) and PSR B 1534-1- 12 jArz oumanian 1995 ; Stairs et al 



2003). iBacker fc SallmenI l|l997l ) claimed erratic emission 



modes from PSR B1821— 2 4 but at only one o bserving 
frequency. In another work, iKramer et al.l l| 19991 ) studied 
PSRs J1022-f-1001 and J1730-2304. In both cases, they de- 
tected profile variations over time-scales of the order of 10 
to 15 minutes; however the data quality for the latter did 
not allow a r igourous statistical analysis. On the other hand, 
iHotan et al.l {2004a) detected no significant variations in the 
pulse profile of PSR J1022-I-1001 and demonstrated that the 
reported profile shape variations could be explained by po- 
larisation calibration errors. 

Small profile changes in PSR J0437— 4715 were de- 
scribed bv IVivo kanan d et ahl l|l998t ) using observations per- 
formed at a very low frequency with only a single po larisa- 
tion. This result was contested bv lSandhu eTall ll 19971). who 
argued that calibration errors were the origin. IVivekanandl 
l|200ll) later argued that the variations are intrinsic to the 
pulsax and correlated with spiky emission in the varying 
component. Variation s in the central regi on of the profile 
were also reported bv lNavarro et all (|l997l ) with 24 minute 
integr ations at 428 MH z but they were not investigated in 
detail. iLiu et al] l|201ll ) developed a sharpness statistic but 
found it insesitive to profile changes in PSR J0437— 4715. 

As described in more detail in ^ this work focuses on 
the stochastic fluctuations in total intensity that arise from 
the subpulse structure observed in single pulses and their ef- 
fect on the timing precision attainable for PSR J0437— 4715. 
In 331 our observations are described along with the applied 
data processing followed by results of timing our observa- 
tions. In ij3]we describe a statistical method useful for de- 
tecting proflle shape variations, then apply it to simulated 
data Eis a demonstration of how it can be used to correct 
ToA residuals. The results of the statistical analysis are pre- 
sented in ^ We summarise our findings and discuss their 
consequences in !j6l which also contains a discussion of other 
possible problems that prevent us from reaching the theo- 
retical timing accuracy. Finally we draw our conclusions in 



2 STOCHASTIC WIDEBAND IMPULSE 
MODULATED SELF-NOISE 

The noise N{t) in equation [1] is normally assumed to be 
dominated by the white radiometer noise. In practice, for 
bright sources and/or high gain antennas, two additional ef- 
fects may contribute significantly. In this section we discuss 
the noise balance for pulsars and demonstrate that the dif- 
ferent contributions are well described by a single statistical 
model. 

Firstly, the noise balance has to include the source itself 
when the flux density of the pulsar approaches the SEFD. 
This noise is intrinsic to the source and is accordingly called 



"self-noise". In the case of PSR J0437— 4715, the mean flux 
at the peak of the pulse profile is of the order of 5 Jy; this 
contributes only ~ 2% to the standard deviation of the to- 
tal intensity, which is domina ted by the SEFD ~ 27 Jy of 
the 20 cm multibeam receiver ^Manchester et cil]|200ll) com- 
monly used for pulsar timing observations at Parkes. Any 
source that can be described as noise (e.g. thermal emission) 
will contribute to the variance of the observed total inten- 
sity of the source. When the signal-to-noise ratio (S/N) is 
low, this contribution is negligible. We note that through- 
out the paper we use S/N values calculated using the noise 
measured in the off-pulse baseline. 

Secondly, dramatic subpulse amplitude modulation is 
a ubiquitous feature of radio pulsar emission (|Rickettlll97"5l ) 
that spans orders of magnitude in i ntensity and duration 
and is a broadband p he nomen on (e.g. Staelin 
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Hankins et al.l 119931 : iHankins fc Eilekl l2007l : IWang et alj 
l2007l ). The s ubpulse emi s sion from PSR J0437-4715 is 
well studied; IJenet et al.l l|l998h observe an exponential 
distribution of peak subpulse intensities with a mean fiux 
density of 16.6 Jy, which is comparable to the SEFD. 
More importantly , the mean subpulse width of 65 ^s 
(|jenet et al.l Il998l ) is about an order of magnitude larger 
than the interval required to sample the mean pulse profile 
of PSR J0437— 4715. Consequently, subpulse intensity 
fiuctuations introduce detectable variations in the average 
pulse profile. Given the stochastic nature of subpulse 
structure, these fluctuations in mean pulse profile shape 
can be expected to introduce significant additional noise in 
derived arrival time estimates. 

In one minute, PSR J0437— 47 15 turns ~ 10^ tim es and 
emits at least ~ 2 x 10^ subpulses (|jenet et al.ll 19981 ). After 
integrating over such a large number of emission events, it is 
no longer practical to consider the impact of an individual 
subpulse. Rather, it becomes necessary to describe the ef- 
fects of subpulse modulation in purely s tatistical term s using 
the fourth moments of the electric field l|Rickettll 19751 ). From 
this perspective, the subpulse modulation phenomenon is 
a noise process tha t contributes to the autocorrelation of 
the total intensity (|Rickettl [l975l ) . For a given source flux 
density, amplitude modulation increases the variance of the 
total intensity and, depending on the time-scale of the mod- 
ulations and the sampling interval of the instrument, intro- 
duces power at non-zero delays in the autocorrelation of the 
total intensity. 

Measured statistical distributions of subpulse intensities 
vary between sources and as a function of pulse longitude 
(e.g. p ower law, log normal, etc.; for an excellent review, see 
ICairns .2004). Regardless of the original distribution, after 
a large number of pulses have been integrated, the central 
limit theorem applies and profile shape variations are well 
described by a multivariate normal distribution. The covari- 
ance matrix that quantifies this distribution contains phase- 
resolved information about the mean autocorrelation of the 
total intensity. 

To summarise, depending on the pulsar's fiux density, 
its emission properties and the used instrument, we can dis- 
tinguish three noise regimes: 

(i) First regime: The pulsar's fiux density is much smaller 
than the SEFD of the instrument. This is the classic regime, 
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in which the noise is temporally uncorrelated between the 
phase bins and homoscedastic (i.e., the variance of noise is 
the same in each phase bin). In this regime the covariance 
matrix of the pulse profiles is well approximated by a diag- 
onal matrix with all the elements on the diagonal equal to 
the variance of SEFD. 

(ii) Second regime: The pulsar's fiux density approaches 
or exceeds the SEFD of the received used. In this regime the 
self-noise cannot be neglected. The noise is still temporally 
uncorrelated between the phase bins, but it is heteroscedas- 
tic; that is the variance of the noise is different in each bin 
and proportional to the sum of squares of the pulsar's flux 
density and the SEFD. The covariance matrix of the data is 
still diagonal, but the non-zero elements are no longer equal. 
In this regime, the on-pulse noise is no longer measured by 
the off-pulse noise and using the latter to calculate the S/N 
can lead to overestimating the achievable timing precision. 

(iii) Third regime: The pulsar is heavily amplitude modu- 
lated with the modulated flux approaching or exceeding the 
SEFD of the instrument. Even though the self-noise con- 
tribution may be negligible, the modulated subpulses can 
approach the SEFD of the receiver, thus contributing signifl- 
cant 'noise' to the averaged pulse proflle. If the sampling rate 
is high enough to resolve the subpulse structure, the noise 
in different phase bins will be heteroscedastic and tempo- 
rally correlated. The off-diagonal elements of the covariance 
matrix will be non-zero in this regime. If the subpulses are 
not resolved, amplitude modulation may still be evident in 
the variation of the modulation index as a function of pulse 
phase. The broadband nature of the impulses will also lead 
to spectral correlation of the noise, which can be detected by 
measuring the covariance of intensity fluctuations in differ- 
ent frequency channels. Therefore, when analysing only the 
covariance matrix SWIMS might be confused with self-noise. 

In this paper, we investigate the effects of the third 
regime, which we call stochastic wideband impulse modu- 
lated self-noise (SWIMS), on pulse arrival time estimation. 
This pulsar-intrinsic noise has also been calle d pulse-phase 
jitter (or jitter noise) , "intermittent emission" l|Gwinn et al.l 
HoTJ) or simply self- noise. We will demonstrate that the 
timing precision of PSR J0437— 4715 is currently limited by 
SWIMS and that its effect on the mean pulse profile is read- 
ily detectable. 

As a function of integration length T, the covariance 
matrix scales as T~^; that is, the effects of SWIMS are re- 
duced by integration, regardless of the dominating source 
of noise. The fact that the covariance matrix scales as 
allows us to study the statistics of single pulses with longer 
integrations. If any pulse-to-pulse correlation were present 
in data, such as arising from drifting subpulses, nulling, 
mode changing, scattering or polarisation calibration errors, 
the scaling of the covariance matrix with integration length 
would deviate from the above proportionality. We note that 
the relative contribution of source-intrinsic noise to the co- 
variance matrix will vary as the flux density of the pulsar 
varies, primarily owing to interstellar scintillation. However, 
after averaging over many scintillation time scales, the rela- 
tive c ontribution of S WIMS compared to the SEFD is con- 
stant (|Kulkarnill 19891 ) : therefore, the relative importance of 
SWIMS is independent of integration length. We note that 
in the first or second regime, noise can be reduced by in- 



creasing the bandwidth. However, because the intensity fluc- 
tuations are typically correlated over wide band widths, the 
noise due to subpulse modulation is not reduced by increas- 
ing the bandwidth and only longer integration times and 
active mitigation can improve timing precision in the third 

regime. 

iGwinn et al.l (|201lf ) recently performed a detailed anal- 
ysis of impulse-modulated self-noise in the context of in- 
terstellar scintillation observations and concluded that self- 
noise may limit pulsar timing precision. In this paper, we 
explore the impact of both temporal and spectral correla- 
tions of intensity fluctuations on pulsar timing and consider 
active mitigation of SWIMS. 



3 OBSERVATIONS AND DATA PROCESSING 

Observations of PSR J0437— 4715 were recorded during one 
week of February 2010 using the Parkes 64 m radio tele- 
scope and the central beam of the 20 cm multibeam re- 
ceiver (|Stavelev-Smith et al.lll996l ). The third generation of 
the Pulsar Digital Filterbank (PDFB3) digitised the volt- 
age data from two orthogonally polarised 256 MHz bands 
and formed 1024 frequency channels using a polyphase fil- 
terbank. After full p olarisation detectio n (following the def- 
initions described bv lvan Straten et ahlboiol ). the data were 
folded at the topocentric period of the pulsar into 1024 phase 
bins. The mean polarisation profile was output every minute 
and a total of 25 hours of data were recorded. The multi- 
beam receiver is equipped with a noise diode that is coupled 
to the receptors and driven with a square wave to inject a 
pulsed polarised reference signal into the feed horn. This sig- 
nal was recorded for three minutes before and after every 64 
minute observation of the pulsar. These data are archived in, 
and can be obtained through, the Australia T elescope On- 
line A rchive and CSIRO Data Access Porta{f| l|Hobbs et al.l 
[20i3). 

The data are stored in the PSRFITS format and all 
processing was done with the PSR CHIVE data processing 
software suite (|Hotan et al ] |2004bl ). First we ensured that a 
recent mo del for the pulsar spi n, astrometric and orbital pa- 
rameters (jVerbiest et al.ll2008l ) was used throughout for our 
data processing. To remove narrow band radio frequency in- 
terference (RFI), median filtering was applied by comparing 
the total flux in each frequency channel with that of its 49 
neighbouring channels. To avoid distortions at the edge of 
the observing band, we rejected five per cent of the frequency 
channels on each side of the band. A search for impulsive RFI 
was performed with the "lawn mower" metho qj. Pulsar ob- 
servat ions were calibrated for polarisation as in Ivan StratenI 
(|2004h . The flux density was calibrated by observing the Hy- 
dra A radio galaxy which is assumed to have a flux of 43.1 Jy 
at 1400 MHz and a spectral index of —0.91. 

A high S/N template for each frequency channel was 
created by integrating the observations obtained during the 
flrst day of data and then used to identify and remove data 
affected by broadband impulsive RFI as follows. First, for 
each frequency channel, the best-fit scale, baseline offset and 

http://datanct.csiro.au/dap/ 
^ http:/ /psrchive. sourceforge.net/manuals/paz 
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Figure 1. The high S/N (15, 000) template for PSR J0437-4715, 
created from 5.5 hours ol observations. The solid black line rep- 
resents the total intensity. 

phase shift (|Tavloj|l993 ) between the profile and template 
were applied to compute the difference between the tem- 
plate and the data. Second, the rms flux of this difference 
was computed after integrating over frequency at the dis- 
persion measure (DM) value of the pulsar and at zero DM. 
RFI will induce a high rms flux at zero DM while, at the 
pulsar's DM, impulsive interference will be smeared across 
multiple phase bins. Hence, if the difference has an rms flux 
value at zero DM that is higher than at the pulsar's DM the 
profile is potentially polluted by RFI. A few hundred profiles 
were examined by eye. The rms ratio at both DMs for the 
RFI polluted difference profiles allowed the determination 
of a threshold ratio above which the remaining profiles were 
automatically tagged as being affected by RFI and rejected 
from further analysis. After all the RFI removal stages we 
were left with 1145 one minute integrations that are consid- 
ered RFI free. From the 5.5 hours of observations taken dur- 
ing the first day of observing, we created a final, frequency 
integrated total intensity template shown in Figure [T] with 
a S/N of 15,000. 

The ToA of each observation not affected by RFI in the 
remaining six days was dete rmined by cross-correlation with 
the template (|Tavloilll993 ). Timing residuals were formed 
from these ToAs and the pulsar m odel using the tempo2 
software package l|Hobbs et al.ll200^ . The ToA residuals for 
these data are shown in Figure [2] as a function of ToA num- 
ber. Note that they are not evenly spaced throughout each 
day. The mean ToA estimation error is only 72 ns and the 
mean S/N is 770. The weighted rms of the timing residu- 
als however is a = 372 ns and the reduced chi squared of 
the fit, /dof , where dof denotes the number of degrees 
of freedom, is 33.8. The unweighted rms timing residual is 
similar: 389 ns. 

The high /dof value could be caused by underesti- 
mation of the ToA uncertainty or because the pulsar model 
does not accurately predict the observed ToAs. To verify 
the estimated arrival time uncertainties, we carried out a 
Monte-Carlo simulation in which each observed profile is 
replaced by an exact copy of the template with a suitable 
amount of white noise added to yield a S/N equal to that of 
each observation when averaged over many realisations. The 
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Figure 3. Comparison of achieved timing precision with the the- 
oretically attainable precision as a function of integration time. 
The line in the dashed region plots the theoretical rms, equal to 
the mean rms obtained from 10^ different realisations of white, 
homoscedastic noise. This values agree with expectations based 
on the width of the mean pulse profile and the S/N . In 95% of 
the simulations, the rms falls within the dashed region. The confi- 
dence interval is much broader at long integration lengths because 
a fixed number of initial pulse profiles is used; hence, at large t, 
fewer independent instances remain, thus biasing the estimate. 



X^ /dof of the timing residuals is always very close to unity 
in these simulations, implying that the ToA uncertainties 
are calculated correctly under the assumption of equation [1] 
In the above simulation, the white noise added to each 
simulated pulse profile is statistically independent of the 
noise in every other profile; in this case, as profiles of roughly 
equal S/N are integrated together, the rms of arrival times 
derived from the integrated totals will be roughly proportial 
to T~^^'^ , where T is the integration length. Consequently, 
the statistical independence of errors in arrival time esti- 
mates is commonly verified using a plot of residual rms as a 
function of integration length, as shown in Figure [3] Here, 
the thick line indicates the mean theoretical expectation 
based on 10^ simulations of 64 ToAs derived from template 
with white noise added. The shaded region shows the 95% 
confidence levels derived from the same simulations. The 
deviation of solid lines from l/\/t behaviour is solely due 
to the small number of points at longer integrations, e.g. 
only 2 points are available at 32 minute integration length. 
The fact that the observed rms follows the expected be- 
haviour suggests that no pulse-to-pulse correlations or anti- 
correlations are important in our dataset. We note that the 
precision measured is far worse than the theoretical expec- 
tation. For example, with 32 minute integrations, we expect 
an rms timing residual of ~ 13 ns but we observe a typical 
value of 52 ns. This factor of ~ 4 worse than the theoretical 
prediction implies that, if this problem was understood and 
fully corrected, then the same observing precision could be 
achieved with integrations 16 times shorter than currently 
required. 

The above simulation does not include any self-noise 
or subpulse modulation; therefore, the predictions in Fig- 
ure[3]are those exp ected from the radiometer equation (e.g. 
iDewev et al.|[l985l l. However, in the real data the variance 
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J0437-4715 (Wrms = 0.372 /zs) post-fit chisq=33.81 




TO A number 

Figure 2. Timing residuals for 6 days of data timed against the standard from Figure[T] The mean ToA estimation error is of the order 
of 72 ns, whereas the weighted rms of the residuals o-j-aj^ is 372 ns. The fit has x/<io/^ of 33.8. For clarity we have plotted the residuals 
as a function of ToA number. The vertical lines are plotted between observations taken on different days. 



of the noise in the off-pulse region understimates the vari- 
ance and completely neglects the temporal correlation of 
the noise in the on-pulse region; that is, the actual noise is 
heteroscedastic and correlated. Increasing the gain of the 
antenna will amplify SWIMS and while the S/N calcu- 
lated using only the off-pulse noise will increase the rms 
timing residual will not decrease. Therefore, longer inte- 
grations will be necessary to achieve a lower rms timing 
residual. More importantly, the SWIMS in the total in- 
tensity from subpulse modulation is spectrally correlated. 
Figure [3] shows a greyscale image of a single pulse from 
PSR J0437— 4715 as function of pulse phase and radio fre- 
quency taken with the ATNF Park es Swinburne Recorder 
(APSR; Ivan Straten fc Bailed l201ll ). The emission clearly 
extends across the entire observed band, producing a high 
degree of spectral correlation of subpulse intensity fluctua- 
tions. To see if To As from independent bands are correlated 
we divided our template and each one minute observation 
into two independent frequency channels and determined 
the timing residuals for each channel separately. Figure [S] 
shows the timing residuals plotted against each other for 
one hour of data processed in this way and shows that the 
ToAs are highly correlated between the two channels (the 
average Pearson product moment correlation coefficient be- 
tween the two sub-bands is 0.91). If the subpulse modulation 
contribution to SWIMS had no impact on the timing resid- 
ual, no such correlation would be present. In addition, the 
rms timing residual of each of the two sub-bands, as well 
as the combined rms timing residual is similar to the value 
obtained for the combined data. Under the assumption that 
the broadband subpulse properties of PSR J0437— 4715 are 
responsible for the scatter in the timing residuals, increas- 
ing the observing bandwidth will not reduce the SWIMS 
component due to subpulse modulation. 

The heteroscedastic and both temporally and spec- 
trally correlated properties of SWIMS can lead to signifi- 
cant statistical bias in arrival time estimates derived from 
sources with amplified flux densities comparable to the sys- 
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Figure 4. Greyscale image of a single pulse from 
PSR J0437— 4715 as a function of pulse phase and frequency 
taken with the APSR instrument. We stress that (a) the sub- 
pulse persists across the whole available band, and (b) that each 
subpulse is very different from the average profile. 



tem equivalent flux density. The following sections report on 
an investigation of one possible method of correcting these 
biases. 



4 METHOD 

As explained in the previous section, failure to account for 
the statistical characteristics of SWIMS leads to measur- 
ment bias and underestimation of arrival time uncertainty. 
In this section, we explore the use of principal component 
analysis (PCA) to correct the statistical bias through a series 
of simplified simulations. These simulations do not model 
the large number of impulsive intensity fiuctuations; rather, 
the simulations demonstrate that the PCA model corrects 
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pulsar timing (see We then form the covariance matrix 
of the dataset by computing the outer product of template 
matched profiles: 



-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 1 
higher frequency band residual [us] 

Figure 5. Timing residuals derived from two independent sub- 
bands every 60 seconds plotted against each otiier to demonstrate 
tiieir remarliable degree of correlation. 



only the arrival time bias due to profile shape variations 
and that no other sources of phase noise ar e incorrectly mit- 
igated . An overview of PCA is given in iHvvarinen et al. 
2001 ) . We extend the analysis introduced by lOemorest 



20071 ) and then present a number of illustrative simulations 



that demonstrate the validity of our method and its im- 
plementation. A very similar methodology has been under 
development by Cordes a nd his collabo rators since the 1990s 
(private communication. ICordes|[l993l V 

The PCA method provides a rigourous and unbiased 
statistical method for analysing temporally correlated vari- 
ations in total intensity. For each one-minute observation of 
PSR J0437-4715 ~ 15 x 10*^ samples are integrated in each 
of the A'^bin pulse phase bins; therefore, by the central limit 
theorem, the fluctuations in total intensity are well described 
by a multivariate normal distribution. If the distribution of 
these fluctuations was strongly non-normal, better perfor- 
mance might be achieved b y a similar method usin g indepen- 
dent component analysis l|Hvvarinen et al]|200ll ). We note 
that the number of pulses integrated in each minute is ap- 
proximately an order of magnitude larger than the number 
of pulses considered in previous studies of profi le stability 
l|Helfand et al.lll975l : iRathnasree fc Rankinlll995^ ■ 

Assume that N such observations have been made of a 
given pulsar. We describe the profile for the i'th observation 
as a column vectoi[f| p'. The j'th element, Pj, is the ampli- 
tude of the j'th bin in the i'th profile. The covariance matrix 
is typically computed after subtracting the mean of all ob- 
servations from each observation. Here, we assume that the 
template, s , is a good estimate of the mean profile and, be- 
fore subtraction, each observation is first adjusted to match 
the template using the best-fit phase shift, scale and offset 
as derived from the template-matching procedure used for 



^ Our notation is defined as follows: All matrices are denoted by 
bold sans serif font (e.g. M). All vectors are denoted by bold italic 
font (e.g. v). An element of a matrix is denoted by Mij. An i'th 
column of a matrix is denoted by M^. An i'th element of a vector 
is denoted by Vi. If there are multiple vectors of a given type, 
we denote the i'th vector by superscripting, e.g. v'. The indices 
i and j always are in range [1,A''] and [1, A'tin], respectively. 



Wi (p" - s) (p" - s) 

1=1 

N 



(2) 



where Wi is the S/N of the i'th profile, the T superscript de- 
notes transposition and the prime superscript signifies that 
the profiles have been matched to the template. The result- 
ing covariance matrix, C, is a symmetric matrix with the 
number of rows and columns equal to the number of bins, 
A^bin, in each profile. We note that at least Abin observations 
are necessary for C to have full rank. Furthermore, the data 
set should be large enough so that all potential modes of 
profile variation are represented. 

Template matching before subtracting the standard 
profile removes three degrees of freedom from the shape 
fiuctuations that are intrinsic to the pulsar signal. For ex- 
ample, to first order, the best-fit phase shift removes all 
variations that correlate with the derivative of the stan- 
dard profile with respect to pulse phase. Removing varia- 
tions with a certain profile shape is equivalent to project- 
ing the Abin-dimensional vector space of the total intensity 
fiuctuations onto the Abin— 1-dimensional subspace that is 
orthogonal to the axis defined by that profile shape. A sig- 
nificant amount of the fiuctuation information may be lost 
by this projection. However, if the best-fit phase shift were 
not first removed, any actual phase shifts would be misin- 
terpreted as shape variations; therefore, this dimension must 
be excluded from the analysis. Similarly, the best-fit scale 
and offset remove variations in pulsar fiux and system tem- 
perature, respectively; these fiuctuations are not the focus of 
this work. In practice, all the data are fit for phase shift, fiux 
scale and baseline offset and these three dimensions are al- 
ways projected out of the available vector space in which the 
pulse profiles are described. The eigenvectors corresponding 
to these three dimensions all have the same eigenvalue (zero) 
and hence together form an eigenspace; we refer to the tem- 
plate matching eigenspace as the fit-space throughout the 
remainder of the paper. 

To characterise the remaining fiuctuations of the total 
intensity, we solve the eigenproblem of the covariance ma- 
trix. The eigenvectors e-* define the principal axes in the 
Abin-dimensional vector space of profile shape variations 
along which the intensity fiuctuations are correlated as a 
function of pulse phase. Sorting the eigenvectors in order of 
decreasing eigenvalue. A, allows us to determine the most 
significant variations. The variance corresponding to each 
eigenvector is equal to the corresponding eigenvalue. 

The eigenvectors form an orthonormal basis onto which 
each residual difference profile can be projected. The pro- 
jection coefficient, Uij, of the i'th difference profile onto j'th 
eigenvector is 



= [P 



(p-s)' 



(3) 



These coefficients can be thought of as the residual of the 
i'th pulse profile in the basis spanned by the eigenvectors e 
and are often referred to as the principal components. 
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After the subspace projection that removes the phase 
shift, scale, and offset dimensions, the remaining projection 
coefficients for each residual profile are uncorrelated; i.e.. 



(4) 



where Sjk is the Kronecker delta. However, these coefficients 
may possibly be correlated with unobservable variations in 
the three dof that have been removed. These correlations are 
exploited by the technique developed in where we in- 
troduce a new method for using these projection coefficients 
to correct the timing residuals for the statistical bias intro- 
duced by pulse shape variability. Simulations to confirm our 
algorithm are presented in H4.2I 

4.1 Correcting the timing residuals: multiple 
regression 

iDemorestI (|2007l ') measured the correlation between the first 
projection coefficient (corresponding to the largest variance 
in the data) and the arrival time residuals. This was sub- 
sequently used to detect corrupted data and he has shown 
that it could be used to remove their deleterious effect on 
the timing residuals. However, his method only used the in- 
formation stored in the first projection coefficient. Here we 
apply a multi-variate statistics method of multiple regres- 
sion to simultaneously remove the effects of multiple varying 
components. 

To predict the statistical bias in ToA estimate, li, 
caused by SWIMS we assume that there is a linear func- 
tion relating this bias to the projection coefficients. We use 
the observed timing residuals, Ri, to determine the best-fit 
parameters for this function. 

We wish to predict li using the linear predictor 



AT 



(5) 



where a and A are the regression coefficients of the linear 
predictor. These are determined from the observed residuals 
by minimising the mean squared error between the predicted 
and observe d residuals: . {Ri — lif ■ The analytic solution 
IS given by (|johnsonlll99a ): 



D 



7 



and 

a — u + A'^ ^ 
where 



D 



{OLi~ < OL >i) {OLi— < OL >i)^ 



iV 



(6) 



(7) 



(8) 



is the covariance matrix of the projection coefficients, 7 is 
the vector of covariances between the residuals and the pro- 
jection coefficients 



7j = 



-<( 



OL ) 

' 7 



> (-R 



(9) 



and 1/ is a vector with each element equal to the mean of 
the observed residuals. The elements [s,i are the mean val- 
ues of OLi and < . > denotes average of vectors. The values 
of A and a allow us to predict the bias in ToA estimation 



induced by pulse shape variations. This predictor has min- 
imum mean squared error and maximum correlation with 
the Ri. Subtracting the predicted bias from the estimated 
arrival times has the potential to reduce the post-fit arrival 
time residual rms. 

The expected improvement in timing residuals can be 
calculated from the projection coefficients and the observed 
residuals as 



where 



P = 



(10) 



(11) 



Here a' is the rms timing residual for ToAs with the bias 
removed and p is called the population multiple correlation 
coefficient. 

ft is necessary to restrict the number of eigenvectors to 
model only pulse variability. Many approaches have been put 
forward in the literatur e for determin ing the number of sig- 
nificant principal axes (|johnsonl |l9"98l). We introduce a new 
parameter, ^j, which is the Pearson's product moment cor- 
relation coefficient between the the timing residuals, i?, and 
the projection coefficients onto the j'th eigenvector, that is 
the j'th row of ol. The standard deviation of non-significant 
^ values is determined. Significant values are identified us- 
ing a Tukey's bi-weighting scheme starting with an initial 
guess of the standard deviation obtained from the median 
absolute deviation. The resulting standard deviation of the 
^ values is a robust and resistant estimator. For more details 
see Andrews (1972) and Hoaglin ct al. (1983). Starting from 
the last correlation coefficient we search for three consecu- 
tive ^ values that are more than three times this measured 
standard deviation. The number of the eigenvectors used 
in all subsequent processing is equal to the index of last of 
these three values, when counting from one. For cases in 
which fewer than three ^ values are significant, we compare 
the results of using only one eigenvector with those of using 
the first five. 

An implementation of this method is publicly available 
as a part of the PSRCHIVE suite. The relevant applica- 
tion is called "psrpca" and it requires the GNU Scientific 
LibrarjH to work. 



4.2 Simulations 

To confirm that our method correctly detects pulse shape 
variations that can be used to correct statistical bias in ToA 
estimation, we carry out three simulations of data with noise 
in the third regime (i.e., SWIMS). In every simulation, each 
observed profile is replaced by a copy of the template profile 
plus white noise and additional varying components. The 
amount of white noise is set such that in many random re- 
alisations of the simulated observation the mean SjN of the 
simulated profile would match that of the given observation. 
Although there can be several subpulses per pulse period, 
we illustrate the technique using a simple model in which 
only a single subpulse is added per minute of observation. 



htt p : / / www . gnu . org/ software /gsl/ 



High S/N 




0.2 0.4 0.6 0.8 1 

phase 



Figure 6. First eigenvector for simulation with white noise, ar- 
bitrary shifts and a von Misses component on top of the template 
profile. Notice that this eigenvector does not look purely like von 
Mises distribution as it has to be orthogonal to the fit-space. 

Note that this procedure does not affect the observing pa- 
rameters (such as frequency, observation time etc.) which are 
held fixed at the values in the actual observations. The re- 
sulting simulated observations are cross-correlated with the 
template to form ToAs (and hence residuals) in exactly the 
same manner as the actual observations. 

We initially tested the trivial cases of simulated data 
with white radiometer noise only with or without arbitrary 
phase shifts applied to the data. As expected, no signifi- 
cant eigenvectors were detected in either case. Attempting 
to correct the residuals regardless of that yields no signif- 
icant improvement in the rms timing residual. These two 
cases demonstrate that our method will not artificially de- 
crease the arrival time residual arising from white noise or 
arbitrary phase shifts. The latter could arise from any phase 
shift such as that due to a gravitational wave and it is im- 
portant that such signal is not removed by the PCA. 

4-2.1 Simulation 1: Single, fixed component 

Pulsar emission is often m odelled as consisting of multiple 
Gaussian components (e.g. iKramer et aLlll999l ). Many pul- 
sars show mode changing where one or more component s are 
active for only a finite amount of time (jWang et al.ll2007l ) . To 
verify that the bias introduced by a single "mode changing" 
component can be detected and corrected, we now include 
an extra component in the profile that varies in amplitude. 
This component is created from a von Mises function (which 
is a periodic analog to a Gaussian distribution) with a nor- 
mally distributed amplitude that has a mean value of zero 
and an rms of 3% of the peak template flux. This compo- 
nent is centred at pulse phase 0.504 and has a concentration 
parameter equal 0.113bins~^. The resulting weighted rms 
timing residual was 255 ns and 'x^ jdoj = 14.4, while with- 
out the additional component the same realisation of white 
noise leads to an rms of 59 ns and a /dof of 1.0. We em- 
phasise that this increased rms residual and high /dof is 
due solely to the pulse shape variations and therefore should 
be detected and corrected using our method. 

We obtain a significant first eigenvalue. The automatic 
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Table 1. Parameters of the multi-component simulation. 

1 2 3 4 5 6 

centre 0.504 0.496 0.512 0.524 0.509 0.520 
amplitude 0.74% 0.36% 0.86% 1.4% 0.7% 0.4% 
K 0.452 0.164 0.098 0.050 0.577 0.104 

determination of useful eigenvectors fails, as there is only 
one significant eigenvector, associated with the introduced 
profile variation. Figure [§] shows the first eigenvector, which 
is different from the introduced component as it is projected 
onto the Afbin — 3 dimensional space, in which the three di- 
mensions corresponding to the fit-space are removed. Al- 
though these degrees of freedom have been removed, the 
residual shape variation is still highly correlated with the 
ToA residual, as discussed above. Correcting the residuals 
using the one significant eigenvector reduces the rms residual 
to 67 ns and the x^ /dof to 1.0. We note that this improve- 
ment agrees very well with prediction from equation [TU] We 
calculate reference residuals by simulating and timing data 
without any additional components but with the same real- 
isation of white noise as in the simulation with the varying 
component. As expected the corrected residuals are highly 
correlated with the reference residuals (correlation coeffi- 
cient of 0.93). This indicates that we have recovered most 
of the signal in the original residuals and have removed the 
effect of the pulse variability. 

We note that after removing the bias in ToA residuals 
the weighted rms timing residual is almost the same as in 
the case of reference residuals; however, the X? /dof can be 
smaller than that of the reference residuals because the ToA 
uncertainties have increased. The total intensity fluctuations 
decrease the cross-correlation between the observation and 
the template, thereby i ncreasing the estimated uncertainty 
(see equation AlO from iTavio^l 19921 ') . 

In order to check the effect of using a different number of 
eigenvectors to correct the timing residuals, we re-analysed 
the data using five eigenvectors. In this case the rms residual 
and x^ /dof were 64ns and 0.92 respectively. The corrected 
residuals are still pleasingly highly correlated with the refer- 
ence residuals. We therefore conclude that our result is not 
highly sensitive to the number of eigenvectors used, but care 
needs to be taken when choosing that number. 

In addition, we also tested an extended version of this 
simulation, where arbitrary phase shifts were included to 
simulate, e.g., the effect of gravitational waves. As expected, 
only the bias in the timing residuals due to the introduced 
proflle variation is removed and the arbitrary phase shifts 
remain unaffected. Again, the residuals are highly correlated 
with reference residuals where the same realisations of white 
noise and arbitrary shifts were introduced. 

4-2.2 Simulation 2: Multiple, fixed components 

The previous simulation dealt with only a single varying 
component while many pulsars can emit radiation simul- 
taneously from several components. Even if only one of the 
multiple components was present in each rotation of the pul- 
sar, several components would be present after integration 
over multiple pulse periods. To demonstrate again that PCA 
and multiple regression do not remove phase shifts that are 
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not caused by pulse shape variations, we now introduce six 
von Mises functions whose amplitudes are allowed to vary 
independently. The parameters of these components are pre- 
sented in Table [TJ where centre is the central phase of a von 
Mises component in pulsar turns; amplitude is the rms of 
the amplitude distribution of given component, in units of 
the template's peak flux; and k is the concentration param- 
eter of the von Mises distribution in units of bin~^. We also 
apply arbitrary phase shifts after adding the varying com- 
ponents. This leads to a weighted rms residual of 372 ns and 
X^/dof of 32.1. Note that the reference residuals obtained 
from a simulation with the same realisation of white noise 
and arbitrary phase shifts, but no additional components, 
have an rms of 206 ns and /dof of 12.2 due to the ar- 
bitrary phase shifts; i.e., we do not expect the rms to be 
of the same order as the ToA measurement error after bias 
removal. 

Our method gives corrected residuals with weighted rms 
residual of 216 ns {x^ /dof — 10.8) and they are highly cor- 
related with the reference residuals. With multiple compo- 
nents the bias is not completely removed because more than 
one projection coefficient correlates with the arrival time 
residual and these projection coefficients may not be statis- 
tically independent of each other. Nevertheless, significant 
improvement in the rms residual and in the X^/dof is ap- 
parent and much of bias is removed. The additional shifts, 
corresponding to unmodelled timing noise processes are still 
unaffected; i.e., the post-correction residuals are highly cor- 
related (correlation coefficient of 0.87) with the reference 
residuals. 

4-2.3 Simulation 3: Single, random component 

We now consider a possibly more realistic case of the pul- 
sar emission being erratic and distributed in phase over the 
whole region in which the average pr ofile is visible; this case 
corresponds to the stochastic nature (jCordes fc Downs! 19851 : 
ICordes|[l993 ) of modulated pulses. In this simulation we al- 
low a single von Mises component per simulated profile to be 
centred anywhere in within the central peak (central phase 
uniformly distributed between 0.479 and 0.518). The concen- 
tration parameter of the component is uniformly distributed 
between 0.055 and 1.386 bin~^. The amplitude of this com- 
ponent is normally distributed with an rms equal to 2.3% 
of the intensity of the template at the centre phase of the 
component. This value is chosen in order that the resulting 
rms timing residual is 380 ns and /dof value of 35.1, both 
very close to the observed values for PSR J0437— 4715. 

Multiple significant eigenvectors are detected, with the 
number varying between different realisations from 10 to 
50. Correcting the residuals reduces the rms to 266 ns and 
X^/do/ to 17.2. We note that in this simulation no arbitrary 
shifts were included and therefore we conclude that the PCA 
method has failed to completely remove the bias in ToAs 
(and hence in residuals) induced by this kind of erratic shape 
variation. As explained in i|4]some fraction of the fluctuation 
power is lost during the vector subspace projection effected 
by removing the best-fit shift, scale, and offset. As before, 
for data integrated over multiple pulse periods, more than 
one varying component per pulse profile would be present 
and this could make it even more difficult to remove the bias 
in ToAs. 
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Figure 7. Distribution of correlation coefficients between the 
residuals and projection onto eigenvectors for the actual obser- 
vations, shown by open circles. Several correlation coefficients are 
more than 3 standard deviations (denoted by the dashed lines) 
above the background noise level. Only the first 150 correlation 
coefficients are shown, the rest of the correlation coefficients look- 
ing very similar as the last shown. The crosses are showing the 
corresponding eingenvalues normalised by their sum, multiplied 
by a factor of 10 and offset by —0.45 for clarity of the figure. 

5 RESULTS 

Applying our method to the observed dataset of 1145 profiles 
leads to: 

• the detection of significant pulse shape variations with 
at least ten significant eigenvectors, 

• a reduction in rms timing residual from 372 ns to 294 ns 
and a reduction in /dof from 33.8 to 21.1. 

The first 150 values of ^ are plotted in Figure [7] which 
indicates that around 10 eigenvectors significantly correlate 
with the observed residuals while the rest are consistent with 
being white noise. The choice of exact number of eigenvec- 
tors would be difficult to make without the rigourous criteria 
described before. In Figure |8] we present the three most sig- 
nificant eigenvectors overlaid on the pulsar template profile. 
The detection of significant eigenvalue-eigenvector pairs is a 
direct consequence of temporally correlated fiuctuations in 
total intensity; i.e., significant shape variations are detected. 
As discussed in [JT] !J5J and 33] it is likely that such varia- 
tions originate from the pulse-to-pulse variations of the pul- 
sar emission. The most significant variation occurs at phase 
0.497, in the main peak of the pulse profile. The majority 
of other significant eigenvectors peak around the main and 
second highest peak in the mean profile. One exception is 
the 10th eigenvector that peaks at phase 0.385, that is in 
the local peak on the left hand side of the main peak. 

Using the 10 most significant eigenvectors to correct the 
bias in arrival time estimates, the rms timing residual was 
reduced by 20% and x^ /dof was reduced by 36%. Only 3% of 
the variance in the timing residuals can be attributed to the 
most significant eigenvector; therefore, multiple regression is 
required to provide the best estimate of statistical bias. The 
third simulation has demonstrated that pulse profile vari- 
ability can introduce bias in ToAs that cannot be removed 
completely using our method and the timing residuals after 
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Figure 8. First three eigenvectors for the observed data. We show only the central part of the profile for clarity. The vertical thick lines 
separate the eigenvectors. The dashed line represents the scaled template added for reference. 



bias correction are still biased. Therefore, even though only 
20% of the rms timing residual is corrected, there need not 
be another explanation for the remaining timing noise. For 
completeness, other effects that might contribute to the rms 
of timing residuals are discussed in the next section. 

To investigate if the intensity fluctuations arise from a 
stationary stochastic process, at least in the wide sense, we 
used five hours of observations of PSR J0437— 4715 during 
July 2009 obtained using the same observing system as de- 
scribed in ij3l These additional data were processed to form 
pulse profiles, arrival times and timing residuals in exactly 
the same way as our main data set. We used the eigenvec- 
tors and regression coefflcients obtained above to correct 
these July 2009 timing residuals. We again achieved a 22% 
decrease in the rms timing residual from 380 ns to 296 ns 
and the x^/dof was reduced by 36% from 39.1 to 24.9. The 
fractional improvement is similar to before, implying that 
the covariance matrix and hence profile variability is sta- 
tionary in time and that after the regression coefflcients and 
eigenvectors have been determined for one dataset they can 
subsequently be applied to other observations of the same 
pulsar obtained with the same instrumentation and observ- 
ing parameters. 



6 DISCUSSION 

We first discuss the method that we have introduced to 
search for and correct pulse shape variability. Second, we 
discuss the astrophysical implications of broad-band profile 
shape variations intrinsic to the pulsar. 



6.1 Discussion of the method 

It is clear from the simulations that our method successfully 
detects significant pulse shape variations and partially cor- 
rects the bias induced in timing residuals due to such vari- 
ations in many cases. The method has certain limitations. 
First, in order for the covariance matrix, C, to have full rank, 
the method requires at least A'^bin observations. With fewer 



observations than phase bins, it may be necessary to average 
adjacent phase bins, or, if prior knowledge suggests that the 
pulse variability occurs in only a restricted region then only 
these bins could be included in the analysis. If full phase 
resolution is required then the covariance matrix can be de- 
termined, but not all o f the eigenvectors can be calculated. 
Alternatively, following iDemores 3 l|2007l ). the PCA method 
can be developed in the frequency domain using only the 
significant harmonics thus reducing the number of required 
observations. 

A large number of observations is desirable for two rea- 
sons. Firstly, dividing a fixed observing time into smaller in- 
tervals provides a greater number of estimates of the pulse 
profile, thereby increasing the S/N of the covariance matrix 
estimate. Secondly, even for A'^ > A^bin our method is lim- 
ited by the SEFD present in all observations. Such noise will 
reduce the precision with which the eigenvectors may be de- 
termined, reducing their ability to fully describe the shape 
variations. For instance, in the timing residuals of simula- 
tions with only white radiometer noise added, even though 
no significant eigenvalues were measured, the rms timing 
residual is reduced by a negligible amount when using just a 
few eigenvectors. When all of the eigenvectors are included 
in the bias removal, the effects of white noise are artificially 
reduced. This occurs because the white radiometer noise af- 
fects the regression coefficients a and A in equation [S] Since 
the eigenvectors are measured using the data that are being 
"corrected', the noise in the eigenvectors correlates with the 
noise in the data. This is similar to the "self-standarding" 
effect that arises when the mean of a set of observed pro- 
files is used as the template to derive arriva l times from 
the sa me data, as described in Appendix A of iHotan et al.l 
(|2005l ). Applying the eigenvectors to a completely indepen- 
dent data set leads to no significant change to the timing 
residuals as the noise in the data no longer correlates with 
noise in the eigenvectors. The degree of correlation between 
individual residual profiles and the white-noise eigenvectors 
may also be reduced by increasing the number of observa- 
tions from which the covariance matrix C is estimated. If 
the eigenvectors are obtained from the data to be corrected. 
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then it is essential to apply the rigourous criteria described 
in i|4.1l to choose only the significant eigenvectors. 

Many pulsar observations are affected by radio fre- 
quency interference (RFI) and the measured eigenvectors are 
extremely sens itive to the pres ence of RFI in the data. As 
pointed out bv lDemoresd (j2007l ). PC A is also a sensitive and 
robust method for detecting RFI and other types of data cor- 
ruption and distortions. It is therefore essential that the data 
used in forming the eigenvectors are unaffected by RFI. In 
our case nearly 32% of the pulse profiles had to be rejected. 
RFI might also be mitigated through the use of a robust es- 
timator of the covariance matri x, such as the minimum co - 
variance determinant estimator (|Hubert fc Debruvnell2ofol ) . 

In contrast to more traditional applications of the PCA 
method, our method relies on first aligning each pulse profile 
to the template using the best-fit phase shift, scale, and base- 
line offset. As a result of this fit, the last three eigenvalues 
are several orders of magnitude smaller (i.e., close to zero) 
and the corresponding eigenvectors spanning the fit-space 
are highly correlated with the template profile, its phase 
derivative, and the baseline offset or their linear combina- 
tions. In other words, the profiles are originally described 
in an A'^bin dimensional vector space and the fit projects the 
data onto A^bin — 3 dimensional vector space in which the fit- 
space has been removed. This removes three degrees of free- 
dom from the remaining eigenvectors and limits the efficacy 
of the correction scheme presented in this paper because 
the fit-space component of the intrinsic shape fiuctuations 
has been removed. Consequently, the bias due to any intrin- 
sic shape fiuctuations that correlate with these eigenvectors 
cannot be fully corrected. This is especially important in 
the case of profile variations correlated with the template 
derivative as such variations will introduce most bias in the 
ToAs. Only the total intensity fiuctuations that are orthog- 
onal to the fit-space contribute to the predictor computed in 
equation [5] The degree to which such fiuctuations are cor- 
rectable depends on how strongly they correlate with these 
three eigenvectors and the degree of correlation between the 
remaining projection coefficients and the arrival time resid- 
ual. 

We note that our work has implications for any relevant 
template matching algorithms. Such algorithms normally as- 
sume that the errors in the measurements of intensity are ho- 
moscedastic and uncorrelated. Our detection of profile vari- 
ations implies that the errors are, in fact, heteroscedastic 
and correlated. We note that including the covariance ma- 
trix, which carries the information about the correlation and 
heteroscedasticity of the noise, into the template matching 
algorithms will yield better estimates of ToA uncertainty. It 
may also remove the necessity of correcting the residuals by 
the means described in this paper because any statistical bi- 
ases caused by the intensity fiuctuations might be removed 
at the time of ToA determination. This is similar to the 
Cholesky decomposition which can remove bias in estima- 
tion of various parameters, such as paral lax, when estimat - 
ing from residuals with red noise present (jColes et al.ll201ll ) . 
Examples of template matching algorithms are th e stan - 
dard ToA derivation algorithm presented bv iTavloil (|l992l ) 
and the matrix template matching algorithm that allows 
all Stokes pa rameters to be used in the ToA estimation 
l|van Straten| [2006'). An unbiased generalisation of template 
matching would be a logical extension of this work. 



6.2 Application to PSR J0437 4715 

Using one-minute integrations of PSR J0437— 4715, we have 
detected shape fluctuations that we attribute to the stochas- 
tic subpulse structure of the pulsar emission. However, for 
many timing applications, the exact cause of the pulse shape 
variations is irrelevant. This technique can be used to cor- 
rect bias and improve sensitivity to any phenomena that do 
not induce shape variations. For instance, in order to place 
a li mit on the existenc e of a gravitational wave background 
(e.g. IJenet et al.ll2006l ) some methods compare the amount 
of power in the timing residuals with the power predicted 
to be induced by gravitational waves. Such waves will not 
affect the pulse profile shape. Therefore, reducing the rms 
timing residuals by accounting for pulse shape variations al- 
lows an improved limit on the existence of a gravitational 
wave signature to be obtained. We note that our correction 
method does not remove any of the signal induced by gravi- 
tational waves or any other phase shift of the pulsar profile, 
as verified by the simulations. 

The long term timing of PSR J0437— 4715 shows signif- 
icant low frequency structure present in the residuals (e.g. 
IVerbiest et al.l |2008| ) . Such red timing noise is a common 
type of non-Gaussianity (in general any asymmetric distri- 
bution will have a similar effect) in the timing residuals. It is 
important to determine the best predictor for correcting the 
residuals based on a short data span, as it will be less affected 
by any non-Gaussian noise. In the presence of non-Gaussian 
components, the best predictor will be affected as the cor- 
relation coefficients ^ between the residuals and projection 
coefficients will be biased toward zero. Presence of Gaussian 
noise (or any other symmetrically distributed noise) will in- 
crease the uncertainty in the estimate of ^ values but will 
not bias them. 

The timing residuals for our Parkes observations are 
only partially corrected by our new method. We have 
demonstrated that partial correction does not necessarily 
imply the existence of other sources of scatter in ToA resid- 
uals that do not affect the pulse shape. At the same time it 
is not possible to completely exclude the existence of other 
effects such as hardware or software errors that are affecting 
the timing residuals at a lower level. Other possible effects 
that can increase the rms timing r e sidual are descri bed in 
detail bv lCordes fc ShannonI l|2010l ): iLiu eTall l|201ll ). Such 
issues are beyond the scope of this paper, which concen- 
trates solely on the pulse shape variability. We note that if 
any other non-white process affecting the residuals could be 
corrected, the residuals induced by profile shape variations 
could be removed more completely as the estimates of ^ are 
affected by the presence of phase noise. 

We stress that the intrinsic shape variations lead to the 
heteroscedastic and correlated component of the SWIM^. 
The uncorrelated component, originating from the standard 
radiometer noise in the weak source limit and the self-noise, 
is described by the diagonal of the covariance matrix C while 
the temporally correlated part is described by the off diag- 
onal elements. However, whatever phenomena contribute to 
the correlated component will naturally also contribute to 

^ We postpone the calculation of the expected amount of ToA 
scatter introduced by SWIMS (based on the measured covariance 
matrix) to following publications. 
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the diagonal of the covariance matrix. As the self-noise and 
subpulse modulation are all measured and described simul- 
taneously by the covariance matrix, it is natural to con- 
sider them all as one phenomenon, which we called SWIMS 
throughout this paper. Two details of covariance matrix we 
would like to stress are that: a) it does not contain any in- 
formation about the spectral correlation of the noise and 
b) the measured covariance matrix also contains a contri- 
bution from the SEFD which can be subtracted if neces- 
sary. Extrinsic sources can also introduce a heteroscedastic 
and correlated component of noise, such as lightning strikes, 
other types of RFI or some instrumental effects, such as 
non-linearity of the backends. The effects of impulsive inter- 
ference should be present across the whole pulse profile as 
they occur randomly in pulse phase. After careful removal 
of RFI, the measured eigenvectors are consistent with white 
noise outside the mean profile (see !j5}. 

As shown in Figure [SI interpretation of the eigen- 
vectors derived from the covariance matrix is complicated 
by the fact that the fit-space has been projected out of the 
A^bin-dimensional space of the shape variations. To investi- 
gate the structure of the shape variations prior to this pro- 
jection, one can make the assumption that the pulsar tim- 
ing model accurately predicts pulse phase (at least over the 
time-scale of the observations) and compute the covariance 
matrix of observed profile residuals after template matching 
by varying only the scale and ofi'set (i.e. no phase shift). In 
this case, the covariance matrix contains the cyclic, phase- 
resolved autocorrelation [autocorrelation function (ACF)] of 
the intensity fluctuations. The mean ACF computed by sum- 
ming elements along the diagonals of this matrix is plotted 
in FigureO The characteristic width of ACF, as determined 
by fitting a Laplace function, is equal to 67 /is, which is 
consistent w ith the aver a ge wi dth of microstructure events 
reported bv Ijenet et al.l ( 19981 ). For comparison, the mean 
ACF formed from the covariance matrix used for bias re- 
moval (i.e. best-flt phase shift removed) is also plotted with 
a dashed line. A large fraction of the fluctuation power is re- 
moved by fltting for the phase shift during template match- 
ing. The phase shift fit removes variations that correlate 
with the template derivative, and the autocorrelation of the 
template derivative reaches its first minimum (below zero) 
at roughly the width of the pulse profile (around 140 /is). 
This corresponds with the first dip seen in the ACF formed 
from data in which no phase shift has been removed; no dip 
is present in the ACF formed after fitting for phase shift. 
The measured ACF may be affected by other non-intrinsic 
effects, especially those related to propag ation through in- 
terstellar medium (ISM: ISmits et aUbo O^. We argue below 
that the ISM is not an important factor in considerations of 
PSR J0437-4715. The ACF shows a periodic ripple (most 
apparent at large phase lags) , which is believed to be an in- 
strumental artefact; it has a period of roughly 100 fis. Its ori- 
gin is currently unknown. Preliminary analysis of data from 
another i nstruinent (Caltech Park es Swinburne Recorder 2, 
CPSR2; iBailedbOOa lHotanll200"^ confirms that this effect 
is intrinsic to PDFB3 as the ACF calculated from CPSR2 
has no periodic ripple present. The pulse profile variations 
are still detected thus confirming their origin as intrinsic to 
the pulsar. 

The interstellar medium can also introduce shape varia- 
tions, such as broadening of the pulse profile, which increases 
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Figure 9. The ACF calculated with the assumption of the timing 
model accurately predicting the pulse phase over the course of 
observations (solid line) and without this assumption, i.e., when 
performing a full template matching (dashed line). The width of 
the ACF is de termined by the characteristic width of single pulses 
(Ric kettlll975l ). The SEFD contributes an unresolved spike at zero 
lag, which is scaled to unit height for the solid line and the dashed 
line uses the same scaling factor. 



quadratically with DM and decreases quartically with fre- 
quency. Given the very low DM — 2.64 cm^'^pc and the 
observing wavelength ~ 20 cm used in this work, the ex- 
pected variations in the pulse width for PSR J0437— 4715 
due to broadening are of the order of 1 ns (|Bhat et al.ll20oj : 
iGwinn et all 120061 ). Interstellar scintillation combined with 
frequency dependence of the pulse profile can also lead to 
fiuctuations of ToAs. As shown in Figure [5]the residuals are 
highly correlated when estimating the ToAs from two sepa- 
rate frequency channels. This high correlation implies that 
diffractive narrow-band interstellar scintillation is not re- 
sponsible for the additional scatter. As we observe variations 
on a time-scale of minutes it is unlikely that broadband re- 
fractive scintillation is a contributing factor to the observed 
pulse shape changes as the time-scale fo r such variations is 
of the order of 1000 s l|Gwinn et al.ll2006h . Based on the high 
degree of correlation between arrival times measured in sepa- 
rate bands (Figure^ and the observed broadband nature of 
single pulses (Figure |4ll, we conclude that the intensity fluc- 
tuations are correlated over wide bandwidths. Consequently, 
increasing the bandwidth does not increase the timing pre- 
cision. Only active mitigation or longer integrations can re- 
duce the timing rms if the fluctuations of total intensity 
are indeed the main cause of the arrival time variations that 
greatly exceed the predicted uncertainty. Simultaneous mul- 
tifrequency observations of PSR J0437— 4715 would help to 
determine if the shape variations are persistent over very 
wide bands. Some observations to date have shown that 
the giant pulses from Crab extend over GHz bandwidths 
(iSaUmen et all 1 19991 : iHankins et"al]|2003l : iHankins fc Eilekl 
200?). Another group has shown that t he single pulses in 
PSR B0329-I-54 persist across 1.3 GHz ( Karas tergiou et al.l 
I2OOII ) but this may not be true for all subpulse structure in 
the general population of pulsars. 

It is also worth considering the impact of polariza- 
tion variations on arrival time estimates. Emission from 
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PSR J0437— 4715 is highly polarised and a sudden change 
in the position angle of the linear polarisation occurs 
near the peak in t he total intensity profile (as noted by 
iNavarro et al.|[r997f) . This implies that poor polarisation cal- 
ib ration will lead to s ignificant profile changes, as quantified 
by Ivan StratenI (|2006l ). For a single dish that is not equatori- 
ally mounted, these variations occur on ti me-scale of hours. 
Pulse profile shape changes detected by ^Vivekana nd et al] 
( 19981) were argued to be caused by calibration errors 
i Sandhuet al]|l997l '). Even though we detect fluctuations 
of the total intensity pulse profile on much shorter time- 
scales, we investigated if this was the case in our data as 
well. We studied the bias in timing residuals induced by the 
measured pulse shape variation as a function of parallac- 
tic angle. We found that there was no dependence between 
these two quantities. We also applied the PCA to uncali- 
brated data and a correlation between the induced bias / 
and the parallactic angle was readily apparent. 

We conclude that the shape variations are more likely to 
originate at the pulsar rather than in the observing hardware 
or from interference. The polarisation calibration has been 
performed sufficiently to alleviate at least minute-time-scale 
fluctuations and does not introduce detectable pulse shape 
variations. The interstellar medium is also unlikely to cause 
such variations. Even if the detected variability is not intrin- 
sic, the presented methodology remains valid. The intrinsic 
variation is expected from the stochastic subpulse structure 
and will be detectable if the pulsar is bright enough. 

Since the detected variations are likely to be intrinsic 
to the pulsar, a question arises whether the profile varia- 
tions in PSR J0437-4715 are related only to SWIMS or if 
they are due to mode changing. We searched for clustering 
in the space spanned by the projection coefficients onto the 
ten significant eigenvectors by applying a friends-of-friends 
algorithm kno wn from n-body si mulations to identify dark 
matter haloes (jPavis et al.lll985h . We did not find any ev- 
idence of clustering in this space and hence conclude that 
the pulse profile variations are not due to mode changing. 

We would like to stress the importance of our work for 
the next generation telescopes, which are likely to provide 
more sensitivity than currently available. With its huge col- 
lecting area of Ikm^, the Square Kilometre Array (SKA) is 
expected to revolutionise pulsar astronomy. One of the Key 
Science Projects of the SKA requires pulsa r observations 
with the highest possible timing precision (|Kramer et al.l 
[iooi). It is assumed that the SKA will observe of the or- 
der of 100 millisecond pulsars with an rms timing precision 
better than 100 ns. With the SKA's phenomenal sensitivity, 
the S/N of a pulse profile should be >1000 on a time-scale of 
only minutes for many pulsars (compared with many hours 
with the Parkes telescope). The short observing times re- 
quired to achieve such high S/N ratios would allow the SKA 
to observe multiple pulsars in a short time. However, the 
increased sensitivity of next-generation telescopes will also 
increase the relative importance of SWIMS as the radiome- 
ter noise is decreased. If the intrinsic pulse shape variations 
that we have detected for PSR J0437— 4715 are typical of 
many MSPs at the observing frequency being used, they 
will induce a floor on timing precision that can be amelio- 
rated only with longer integration times and active mitiga- 
tion using methods suc h as the one presented in this work. 
ICordes fc Downs! (|l985h demonstrated that, for the major- 



ity of their sample of 24 pulsars, timing precision is likely 
to be limited by pheise jitter. No millisecond pul sars were 
included in this sample. Ishannon &: CordesI l|2010l ') later ar- 
gue that the timing noise of millisecond pulsars is similar 
to that of classical pulsars, only much smaller. Although 
SWIMS has not yet been detected in most pulsars, it is likely 
to be revealed with better instrumentation, more sensitive 
telescopes or longer data spans. Since the detected pulse 
shape variations are likely to be very broadband, increasing 
bandwidth will not reduce the bias introduced by SWIMS. 
This must be considered when predicting the potential sci- 
ence of current and future pulsar timing array projects and 
the observing time and strategy necessary to achieve the 
stated goals. For example, with an SKA-like telescope, if 
many pulsars in a timing campaign are limited by SWIMS, 
then it is better to observe multiple pulsars simultaneously 
with fraction of the array for longer time rather than us- 
ing full sensitivity to observe pulsars one by one for a short 
time. Some proposed astrophysical experiments demand ex- 
tremely accurate ToAs over short intervals of the order of 
minutes, such as when a pulsar passes behind black hole in 
a close binary. SWIMS will make such experiments difficult 
or impossible. 

The relative importance of the correlated component 
of SWIMS is also expected to vary between pulsars as it 
depends on intrinsic subpulse emission properties and the 
shape of the average pulse profile. The fractional improve- 
ment in rms timing residual is expected to vary from case 
to case and it can be hoped that for pulsars with simpler 
and/or narrower profiles, variability in the subpulse struc- 
ture will be less severe. As demonstrated by our first simu- 
lation, in simple cases our method works very well and can 
completely remove the statistical bias in ToAs for some pul- 
sars. Our method can be used to identify pulse profile modes 
which can lead to improved timing. 

We note that, for current telescopes, equation [T] is a 
good approximation for the vast majority of MSPs, which 
have time averaged mean flux densities an order of mag- 
nitude smaller than PSR J0437— 4715. For example, in the 
Parkes Pulsar Timing Array sample, fluxes v ary between 1.3 
and 13.8 mjy (see Table 2 of lYan et al.ll201ll ) with a median 
value of 2.4 mJy, compared to the mean value of 150 mJy for 
PSR J0437— 4715. Consequently, the profile variability aris- 
ing from SWIMS has been neglected to date. We note that 
future telescopes like the SKA are likely to perform timing 
array experiments at higher frequencies to avoid some of 
the problems caused by the interstellar medium. Whether 
SWIMS will be a crucial limitation to precision timing for 
the majority of pulsars at all observing frequencies remains 
to be seen. 



7 CONCLUSION 

We have developed an extended principal component anal- 
ysis method that is applicable to searching for pulse 
shape variations in pulse profiles. Applying this method to 
PSR J0437— 4715 shows the presence of pulse profile vari- 
ability that is likely to be intrinsic to the pulsar. The statis- 
tics of this variability are consistent over many months. The 
detection of significant intensity fluctuations implies that 
self-noise may be a limiting factor for timing precision of 
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PSR J0437— 4715 for current generation of telescopes. Fu- 
ture technological developments including construction of 
larger antenna and increased instrumental bandwidth will 
not improve timing precision as the subpulse structure is 
a source-intrinsic broadband phenomenon. However, the ef- 
fects of SWIMS can be partially corrected by the method 
presented in this work and the proposed generalised tem- 
plate matching. 
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